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Abstract 


This  paper  addresses  orbital  perturbations  in  the 
two-body  problem  of  an  earth-orbi ting  solar  sail  undergoing 
free  coning  motion  in  a  circular  orbit.  The  coning  motion 
controls  both  the  magnitude  and  direction  of  the  solar 
radiation  pressure  force.  The  equations  of  motion  are 
expanded  from  the  Lagrange  planetary  equations  in  their 
acceleration  component  form  and  are  solved  assuming  small 
changes  in  the  orbital  elements  over  the  period  of  one 
orbit,  n 

Resonances  are  observed  between  the  mean  motion  of  the 
sail  and  its  precession  rate.  The  one-to-one  resonance  case 
is  examined  and  numerical  methods  are  employed  to  verify  the 
analytic  results  for  the  circular  orbit  case. 
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ORBITAL  MOTION  OF  A  FREELY.  CONING  SOLAR  SAIL 


I.  Introduction  and  Problam  Statement 


Introduction 


Spacecraft  propall  ad  by  solar  radiation  prassure  ara 
baing  considerad  -for  usa  in  spaca  exploration.  Sucl 
vahi  cl  as  raquira  large  areas  of  axtramaly  light-weight, 
highly  ref lectiva  sails  to  convert  the  linear  momentum  o\ 
sunlight  into  thrust  for  tha  spacecraft.  While  the 
acceleration  of  such  a  vehicle  is  quite  small,  it  is 
continuous,  variable,  controllable,  and  does  not  require  the 
expulsion  of  mass  overboard  as  with  traditional  momentum 
exchange  techniques. 


Solar  sailing  vehicles  can  provide  an  adequate  means  of 
transportation  across  the  long  distances  typically  found  in 
space  travel.  The  force  applied  to  feasible  sails  is  large 
compared  with  space  drag  (Ref  8).  While  the  accelerations  of 
such  vehicles  are  small,  the  velocity  built  up  over  time  is 
substantial.  A  solar  sailing  vehicle  which  can  achieve  a 
velocity  sufficient  to  escape  the  solar  system  is  feasible 
(Ref  4).  Escape  from  Earth  orbit  is  possible  after  a  month 
or  two  (Ref  4).  The  solar  sail  has  the  additional  advantage 
of  simplicity  of  design  when  compared  with  other  low-thrust 
propulsion  systems  (Ref  8). 


In  attempts  to  determine  the  best  steering  program  tor 


solar  sails,  a  variety  of  approaches  have  been  tried. 
Garwin  (Ref  4)  began  the  investigations  with  a  sail  that 
could  be  furled  when  sailing  into  the  Sun’s  rays,  and 
unfurled  when  sailing  with  the  solar  flux.  This  is  similar 
to  a  method  proposed  by  Van  der  Ha  and  Modi  (Ref  9>  which 
involved  "on— off  switching".  London  (Ref  6)  sought  the  one 
best  constant  angle  for  a  given  mission.  A  more  advanced 
method  was  introduced  by  Fimple  (Ref  3)  who  determined  a 
steering  program  based  on  the  time  rate  of  change  of  energy 
increase.  All  of  these  approaches  have  in  common  that  if 
the  sail  is  rotated  at  all  it  is  rotated  with  the  angular 
momentum  vector  coinciding  with  a  body  axis  of  symmetry, 
i.e. ,  spinning.  This  article  represents  a  departure  from 
this  traditional  approach. 

The  question  is,  essentially,  how  to  best  control  the 
thrust  of  a  solar  sailing  vehicle.  Because  of  the  inherent 
flimsiness  of  light-weight  sails,  quick  maneuvers  are  not 
possible  with  solar  sailing  vehicles.  Looked  at  another 
way,  if  a  method  can  be  found  which  both  provides  adeauate 
performance  and  restricts  the  sail  to  smooth  rotations,  then 
the  necessary  structural  strength  and  weight  of  the  sails 
can  be  minimized.  Adequate  performance  is  defined  in  terms 
of  the  spacecraft’s  ability  to  increase  the  semimajor  axis 
of  its  orbit  and  to  change  the  inclination  of  its  orbit. 
Specific  performance  criteria  depend  on  mission 


requirements.  Control  o-f  the  coning  motion  would  be 
accomplished  by  control  o-f  the  magnitude  of  the  angular 
momentum  vector. 

This  work  uses  techniques  pioneered  in  Capt  Salvatore 
Alfano’s  thesis:  "Low  Thrust  Orbit  Transfer"  (Ref  1).  For 
those  not  familiar  with  the  work,  the  consideration  of 
orbital  changes  for  very  low  thrust  vehicles  is  broken  into 
two  parts.  The  first  part  is  a  fast  timescale  problem  in 
which  many  orbital  parameters  can  be  considered  constant 
over  one  orbit,  the  perturbations  being  found  and  added  at 
the  end  of  the  orbit.  The  second  part  is  a  slow  timescale 
problem  in  which  many  perturbed  orbits  are  linked  to  achieve 
a  transfer.  What  follows  here  is  the  derivation  and 
examination  of  the  equations  for  the  fast  timescale  problem 
for  solar  sails. 

Problem  Statement 

The  objective  is  to  derive  the  perturbation  equations 
for  an  earth-orbiting  solar  sail  which  is  undergoing  a 
coning  motion.  The  sail  is  modeled  as  a  flat,  rigid,  plate 
which  is  undergoing  the  free  coning  motion  of  an 
axisymmetric  rigid  body.  During  the  course  of  one  orbit, 
the  sail  is  considered  to  be  torque  free.  The  earth  is 
modeled  as  a  point  mass.  Because  of  the  low  thrust  of  the 
vehicle,  only  small  changes  in  the  orbit  are  possible. 


Thus,  the  semimajor  axis  <a> ,  inclination  (i),  eccentricitv 
(e)  ,  argument  a-f  periapsis  <6J)  ,  and  the  longitude  o i  the 
ascending  mode  (A)  are  considered  constant  over  the  orbit 
and  then  updated. 

The  perturbation  equations  are  to  be  examined 
analytically  to  determine  what,  i-f  any,  important  resonances 
arise.  A  promising  resonance  case  will  be  examined. 


II.  Attitude  Dynamics  of  a  Freely  Coning  Solar  Sail 


Fundamentals _of  Coning  Motion 

The  attitude  motion  of  the  sail  is  that  o-f  the  freely 
coning  motion  of  an  ax i symmetric  body.  The  body  is  a  flat, 
rigid,  plate  which  exhibits  perfectly  specular  reflection  on 
both  sides.  A  review  of  the  fundamentals  of  free  coning 
motion  is  appropriate. 


Coning 

motion 

is 

the 

motion  of 

a 

torque  free 

ax i symmetric 

rigid 

body 

in 

the  case 

where 

the 

angular 

velocity  vector,  cu 

,  is  not 

col inear 

with 

the 

angular 

A 

momentum  vector,  H.  The  assumption  that  the  body  is  torque 

A 

free  implies  that  H  is  constant  in  direction  and  magnitude. 

A  A 

It  has  been  shown  (Ref  S)  that  H  is  coplanar  with  to  and  a 

A 

vector  along  the  body's  axis  of  symmetry,  b3  .  Furthermore, 
since  the  orientation  of  the  angular  momentum  vector  is 
constant  in  inertial  space,  the  plane  which  contains  the 
three  vectors  rotates  about  the  angular  momentum  vector  (Fig 

A  A 

1).  We  can  designate  the  angle  between  M  and  b3  as  the 
nutation  angle,  0 ,  which  is  constant.  A  second  angle,  rp>  , 

A  A 

is  the  angle  between  CJ  and  b3  .  This  plane  rotates  with  a 
constant  angular  velocity  A.  As  the  vector  plane  rotates, 

A 

describes  a  space  cone  around  H.  This  cone  remains  fixed  in 

A  A 

inertial  space  and  a  cone  defined  by  b3  and  to  rolls  without 
slipping  on  the  space  cone.  For  a  body  with  moments  of 

-  5  - 


inertia  A’,  A*  and  C’ ,  where  C’  is  the  moment  of  inertia 

A 

about  b3  ,  the  -following  relationship  applies  (Ref  5): 


tan  6*  =  Atan^ 

C' 

^  (i) 

For  cases  where  C'  >  A* ,  this  implies  that  as  shown 

in  Figure  1. 

A 

Consider  that  H  is  col inear  with  an  inertial  reference 

A  A  A 

axis,  K,  and  that  the  I  and  J  axes  -form  a  right-handed 

A 

orthogonal  set  with  K  (Fig  2).  Consider  also  a  projection  o-f 

A  A  A  A 

b  onto  the  I,J  plane.  This  projection  will  rotate  about  K 

3 

A 

and,  therefore,  about  H  at  the  constant  precession  rate  V 
given  by* 


~y  ■  Cot /(A- C')  cos 0  <2) 

A 

where  oC  is  the  angular  velocity  of  the  body  about  the  b3 
axis  (Ref  5).  For  1/  to  be  positive  in  the  case  where  C’  > 
A’ ,  oC  must  be  negative.  Since  (X  is  just  the  spin  of  the 
plate  about  its  axis  of  symmetry,  this  presents  no 
difficulty. 


The  main  points  of  interest  are  that  the  precession 
rate,  nutation  angle,  and  the  angular  momentum  vector  are 
all  constant.  These  facts  provide  the  basis  for 


determination  of  the  magnitude  and  direction  of  the  thrust 


vector.  From  this  general  groundwork,  specific  equations 
for  the  attitude  motion  of  the  solar  sail  can  be  found. 


Attitude  Motion  of  the  Sail 

The  solar  sail  is  a  flat  plate,  perfectly  specularly 
reflective  on  both  sides  and  with  moments  of  inertia  A’,  A’, 


and  C* .  C’ 

is 

the 

major 

axis  of 

inertia  and  lies 

perpendi cul ar 

to 

the 

plane 

A 

of  the  sai 1 . 

A  A 

The  body  axes  of 

A  A 

the  sail  are  designated  b3  ,  ,  and  b,  ,  where  b(  and  bz 

A 

are  axes  about  which  the  moments  of  inertia  are  A’,  and  b3 
is  the  axis  about  which  the  moment  of  inertia  is  C’ .  (Sec 
Fig  3.) 

The  original  orbit  plane  will  be  considered  inertial. 

A  A  A  A 

The  inertial  axes  will  be  designated  i,  j ,  and  k  where  k  is 

A 

normal  to  the  orbit  plane,  i  is  in  the  orbit  plane  and  in  a 
plane  perpendicular  to  the  ecliptic  and  containing  the 
Earth-Sun  line,  and  j  forms  a  right-handed  orthogonal  set  as 
in  Fig  4.  The  angle  between  the  first  point  of  Aries  and  the 
Earth-Sun  line  is  considered  to  be  constant  over  one  orbit. 

Both  the  magnitude  and  direction  of  the  thrust  are 

A 

functions  of  time.  The  magnitude  in  the  -b  direction  is 
simply: 


T  =  5hAg(b3-s)2 


<3) 


where  s  is  the  solar  constant,  Pi  is  a  sail  parameter 
dependent  upon  the  re-f lecti vi ty  and  transmissivity  o-f  the 

sail,  A  is  the  area  o-f  the  sail,  A(b  •  s)  is  the  projected 

area  o-f  the  sail  on  a  plane  perpendicular  to  the  Sun’s  rays, 
and  g  is  the  unit  conversion  factor.  The  unit  vector  s, 

shown  in  Fig  5,  gives  the  direction  of  the  Sun  and  is 

defined  by 


/'  •  V 

S  =  COS  L  L 


A 

sinl  k 


(4) 


The  definition  of  i  used  here  and  shown  in  Fig  4 
implies  that  XI  is  fixed  at  a  negative  one-half  pi  radians 

A 

from  the  i  axis.  There  is  great  simplification  in  this 
choice,  and  little  loss  of  generality  for  the  current 
problem.  The  maximum  change  in  inclination  can  still  be 

seen. 

A  A 

The  direction  of  the  thrust  is  -b}  ,  assuming  that  b^ 

A  A 

is  on  the  sunward  side  of  the  sail.  Since  i,  i,  k,  s,  n  and 
A  are  known  and  constant  over  one  orbit  it  remains  only  to 

A 

find  b3  expressed  in  terms  of  functions  of  time. 

A 

Given  that  the  angular  momentum  vector  H  has  an 


arbitrary  orientation  in  inertial  space,  we  can  define  a 


a' 


a* 

v  • % 


r> 

i 


m 

>.V.' 

•  •  fc  ' 


,-.  J* 

i 

>.V. 

$ 

Xv 

m 


C<: 

■:■!■. 


AAA  A  A  a 

coordinate  system  I,  J,  K  such  that  K  is  along  H,  I  is 

A  A  a  A 

perpendicular  to  H  and  lies  in  the  i ,  j  plane,  and  J  -forms  a 

right-handed  orthogonal  set.  This  permits  the  definition  of 
a  . 

H  by  means  of  two  angles  of  rotation  j  and  7  ,  as  shown  in 

Figure  6  (Ref  10).  Any  vector  in  this  I,  J,  K  frame  can  be 

A  a  A 

expressed  in  the  i,  j,  k  inertial  frame  by  the 
transf or mat ion 


where 


R.(x) 


r. 

L 

I 

' 

j 

>  = 

R3(-f)R,C 

-r})< 

J 

> 

k 

K 

k  J 

f'  0 

0 

cosx  sinx  ol 

0  cosx 

sinx 

R3(x) 

-sinx 

cosx  0 

0  -sinx 

COSX 

— 

0 

0  1 

(5) 


(6) 


The  unit  vector 


A 

b 


AAA 

is  defined  in  the  I,  J,  K  frame 


shown  in  Fig  7.  This  leads  to  an  expression  for  b3 
form* 


in  the 


A  A  A  A 

b3=  sin@cos(v(t-t0))I  +  5in^sin(-^(t-te))  T+ 


cosi9K  <7) 
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The  nutation  angle  0  is  a  constant  over  one  orbit.  The 
precession  angle  "V<t-t0  )  is  obviously  a  function  of  time. 


For  simplicity,  t0  is  set  at  zero  and  a  phase  angle,  f" ,  is 

A 

used  to  offset  the  starting  point  away  from  the  I  axis  when 
desirable.  As  will  be  seen,  phase  angles  are  crucial 
because  a  given  coning  motion  can  create  different  effects 
depending  on  initial  conditions. 

A  f\ 

To  get  thrust  magnitude  we  need  (b^*  s)  which  leads  to 

A  A  A  A 

the  necessity  of  expressing  b  in  terms  of  the  i <  j,  k 

o 

reference.  This  is  accomplished  by  transf ormati on  Eq  5, 
resulting  in 


L  J 

r  > 

L 

f 

sOcZ'tc.^  -  sesistcysf  +c9sr}tf 

j 

. 

■ 

+S^svtc/?c^  -  c > 

k 

^  / 

S  9svtsij  +  cd  cq 

(8) 


which  implies  that 


A 


A 

S 


d,cos^f+^)- d^sinC^t+T)  +  d3 


(9) 


where  constants  d(  ,  d„  ,  and  d3  are  defined  by: 


d(  =  sin  9cos%cosi 


d2  =  sm G  cosr|sia^ cosi  +  sin^sinrjsinl 
d  =  cos#  sini^sin)?  cost  -  cos#cosrjsml 


It  is  convenient  to  define  a  set  of  axes  with  their 
origin  at  the  center  of  mass  of  the  spacecraft  and  with  axes 

A  A  A 

U,  V,  and  W  in  the  orbit  radial,  tangential,  and  normal 
directions,  respectively.  To  get  the  direction  of  the 

A  AAA 

thrust,  b  must  be  expressed  in  the  U,  V,  W  reference 
frame.  This  is  accomplished  by  noting  Figure  8  and  writing 
the  transf ormation 


v  =  R3(f)  j  =  R3(4)R3(-^)RI('-f?)  J 

w)  l  kJ  Ik, 


where  f  is  the  true  anomaly  and  ^  is  a  phase  angle. 


Transf ormation  <11)  results  in 


U  flcvd  -  /4svcJ  + /33c$  + +  /3S  sz/sf  +/|s^ 
b3  <  V  >  =  <  -^3,  cvsf  +/^svsf  -^3  5^  +/5yC^cf  +^5  ST^C-5- +^cf 

W  /?7  57/  +  4 

\  J  V  , 

sssin,  cos  ,  7/=i/t+Y,  fsnt  +  ^ 

where  the  constants  A  through  A  are  defined  by: 


=  sin^cosJj 
^  =  sin0cosi|sin  J 
=  cos©  sinrjsin^ 
y6^  =  si  n<9  smf 
»  sm0  cosr| 

*  coS©simjcos£ 

i S-j  =  sin©  sinr^ 

^  *  COS&COSlJ 


These  results  provide  the  accelerations  in 
tangential,  and  orbit  normal  directions. 
Newton’s  Second  Law: 


F=  ma  =$  a=  F/m 


then 


(13) 


the  radial, 
Remembering 


(14) 


U  =T/m  =  D(b3-s)l(t.U) 


III.  The  Orbital  Equations  of  Motion 


The  objective  is  to  develop  the  orbital  equations  of 
motion  from  the  Lagrange  planetary  equations  in  their 
acceleration  component  form.  The  Lagrangian  equations  are: 


X.Uesin-5-  +  ZV  ( \  +•  ecos£ ) 

at 

n(l  -e^V4 

di  • 

W  cos  (Kw)  ( l  - 

■*)* 

at 

na(  l  +  ecosf , 

) 

cLo.  = 

W  sin (f  +  uj)  ( i- 

dt 

na(  1+  ecos$) 

sine 

do)  - 

-UcosHl-e1)'"  +  V(l-rfl(^ecw^ 

0  sin f 

at 

nae 

nae  ( i  +  e  cos-^ ) 

-Wsmtf  +  aj)  cote  f  l- 

no.  (  i  + 

e  cos$ ) 

i§  = 

U_0-  e*V'*sinf 

dt 

na. 

r  i 

+  V  (l-e1)1'1 

1 1  ecos$  -  i  -  el 

nae 

1+  ecosf 

(19) 


(20) 


(21) 


(22) 


(23) 


where  U,  V,  and  W  are  the  spacecraft  accelerations  in  the 
radial,  tangential,  and  orbit  normal  directions  and  n  is  the 
mean  motion.  The  true  anomaly  is  designated  f  and  consists 


■f  -  n(t-t0)  +0 


where  t0  will  be  taken  to  be  zero  and  ^  is  a  phase  angle 

A 

between  the  i  axis  and  the  beginning  of  the  orbit. 


For  circular  orbits,  cj  is  arbitrarily  set  to  zero,  and 
the  time  derivative  of  the  argument  of  periapsis  is 
undefined.  This  reduces  the  number  of  equations  of  motion 
to  four.  Simplifying  for  e  equals  zero  we  gets 


cia  =  2V 
dt  n 


de  =  U  sirrf  +  2Vcos~f 
dt  na  na 


di  =  W  cos~f 
dt  na 


L* 
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These  are  the  orbital  equations  of  motion  -for  the 


freely  coning  solar  sail  in  the  Lagrangian  acceleration 
component  form. 

To  obtain  the  changes  in  orbital  parameters,  Eqs  29-32 
are  integrated  with  respect  to  time  over  one  orbit,  the 
period  being  determined  from 


(33) 


Perturbation  Solutions 

To  obtain  the  general  perturbation  solution  for  the 
orbit  of  a  freely  coning  solar  sail,  Eqs  9  and  12  are 
substituted  into  Eqs  29-32  which  are  integrated  with  respect 
to  time  between  the  limits  zero  and  Ti  Since  these  equations 
contain  both  f  and  ~v  which  are  functions  of  time,  the 
perturbations  solutions  are  quite  lengthy.  For  this  reason 
they  have  been  placed  in  Appendix  A.  The  reader  is  strongly 
urged  to  examine  them  in  some  detail  before  proceeding.  The 
perturbation  equations  are  given  here  in  abbreviated  form 
for  notational  clarity: 


.T 


Aq  =  \  da  dt 

)  dt 


(34) 


IV.  Resonances 


Inspection  of  the  equations  in  Appendix  A  indicate  the 
presence  of  resonances.  These  appear  as  terms  in  the 
denominators  which  could  go  to  zero  if  the  correct 
relationship  between  "V  and  n  existed.  When  a  resonance 
condition  is  met,  the  term  in  which  the  resonance  occurs 
becomes  large.  This  is  the  key  to  maximizing  the  small 
changes  in  the  orbital  elements,  which  is  our  objective. 

Several  different  resonances  occur  in  the  perturbation 
solutions.  In  order  to  more  easily  select  useful 
resonances,  all  resonances  have  been  displayed  in  Tables 
I— IV.  The  resonances  are  given  in  the  form: 

(38) 


where  x  is  the  resonance  number  given  in  the  table.  The 
tables  display  the  resonances  in  a  grid  which  indicates 
which  constants  each  resonance  affects.  For  instance,  the 
one-to-one  resonance  enlarges  the  constants^,  and 
^5$  each  multiplied  by  the  constant  dg ,  in  the  A  a 
equation.  A  negative  sign  beside  a  resonance  number  means 
that  this  term  tends  to  reduce  the  perturbation,  assuming 
both  phase  angles  ( ?  &  Y'i  are  zero.  Adjustment  of  the  phase 
angles  can  change  the  sign  of  a  given  term.  Using  these 


TABLE  II 


Resonance  tar ms  in  a 


(Constants)  d*  2d,  da  d^ 


A 

-2,  -2/3 

2,  2/3 

-2,  2/3 

A 

2,  2/3 

-2,  2/3 

2,  -2/3 

A> 

-1 

1 

1 

At 

2,  2/3 

2,  2/3 

2,  -2/3 

As 

-2,  -2/3 

-2,  2/3 

-2,  2/3 

A 

1  * 

1 

-1  * 

(Constants) 


2d,  d 


2d,d3 


TABLE  III 


Resonance  terms  in  i 


(Constants) 

< 

2d  d 

i  i 

4 

A7 

-1,  -1/3 

-1,  1/3 

-1,  1/3 

A% 

1/2 

1/2 

-1/2 

(Constants) 

2d,d3 

2dld3 

<4 

A 

1/2 

1/2 

-1 

A 

1  1 

TABLE  IV 

Resonance  terms  in 

0 

(Constants) 

2d  d 

\  a. 

4 

; 

A 

1,  1/3 

-1,  1/3 

1,  -1/3 

A* 

1/2 

1/2 

-1/2 

r 


(Constants) 


2d.  d 


2d  d 


tables,  a  useful  resonance  can  be  chosen. 


Qne7tO30ne  Resonance  Case 

The  perturbation  solutions  for  the  one-to-one  resonance 
case,  integrated  over  one  orbital  period,  are  given  here. 

A  a  =  a3/2DT,{(d,/tt)(  {.//31+3/3l)cos(0-y)  -  (3/5l+/55)sin(^-‘y')} 

+  (//4>s-yd’,)  cosip-V')}  +  Zd,d3{ 

/^co s((f>-Y)  s in(fl-Y')}  +  Zc\^3i /S^cosLefi-Y)  si  n(^ -¥'.)] 

+(dz/i){(3ySL+/^’|)cos(^-'^)  -  (/3^  3/35)  s\n{(p~Y)}  + 

d3  {(/S^  +/3i)cos  Lp-V)  -  ( yf,-*-/^)  sinlfi-V')} 3 

(39) 

Ae-(a*D7/8)[<!*[y%  5in(z(p-y'))  6  +  cos  (zc^-?)))} 

+^Ja  cos 

+A(  6  -  cos  (zC^-VO))]  +  ZdId3{(y5z+/54)cos(z(^-y'))  -  (//3l  + 

/S5  )sin(z(^-y>))  +  6/5,}  +  Zdzd3{(  ydi+jdf)  sin(z(ff-r))  +  (/$, 

+/55)  cos  (l.(fi-V'))  -  6y65}  +  !Zd3/56} 

(40) 


Al  =  (a^DT/z){2d3/58{d1  cos(^-V')  +  djrSinC^-y-)}  " 

/?,  {[(df/‘t)  +  f3djA)  +  dj]  sin^-VO  +  (d,dz/z)  cos  (*-*>}} 

(41) 
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Ail  =  (  a/zDTAsini){ ^  {  [  (  df/l)  +  +  cosip-f) 


—  (d, <dz/2.)  sin  (^— *2-/3 s  cj3 {  d t  sin  (fl-Y)  -  ciscos  (.$ -?■)}} 

(42) 

The  one-to-one  resonance  case  has  several  interesting 
properties.  First,  it  the  most  -frequently  encountered 
resonance  case  in  the  A  a,  Ai»  and  AH  equations.  Since 
the  performance  requirement  is  to  have  a  large  change  in  the 
semimajar  axis  and/or  the  inclination,  the  one-to-one 
resonance  case  is  pre-f erred.  Secondly,  when  the  one-to-one 
resonance  is  selected  and  the  integration  is  taken  over  one 
orbital  period,  i.e., 

nT  =  Zl T  (43) 

then  virtually  all  terms  not  containing  a  one-to-one 
resonance  disappear.  The  only  exception  is  the  last  term  in 
the  A  e  equation  which,  along  with  four  other  Ae  terms 
which  are  one-to-one  resonant,  cannot  be  integrated  to  zero 
under  any  circumstances.  These  five  terms  are  non-periodic 
and  can  be  eliminated  only  by  proper  selection  of  constant 
angles.  Non-periodic  terms  occur  only  in  the  Ae  equation 


and  are  indicated  by  an  asterisk  in  Table  II. 


Thirdly,  the  one-to-one  resonance  case  is  the  only 
integer  resonance  case  in  the  Aa,  A  i  r  and  A£2  equations. 
To  integrate  non-resonance  terms  to  zero  -for  the  -fractional 
resonances  requires  a  period  o-f  more  than  one  orbit.  This 
limits  the  orbits  -for  which  the  perturbation  assumptions  are 
valid  to  smaller  semi major  axes  than  for  the  one-to-one 
resonance  case.  Conversely,  if  one  orbit  is  taken  as  the 
upper  limit  of  integration,  none  of  the  terms  would  be 
eliminated  and  the  simplicity  of  the  analytic  method  would 
be  destroyed. 


Other  Resonances 

The  two-to-one  resonance  case  changes  e  only.  In  the 
event  that  the  mission  objective  involves  changing  the 
eccentricity  only,  this  would  be  the  resonance  of  choice. 
Because  of  the  circular  orbit  assumptions,  this  case  is  not 
investigated  further  here. 

The  one-half -to-one  resonance  case  would  be  the 
secondary  choice  for  changing  a  and  i,  but  it  should  be 
noted  that  the  two-to-one  resonance  terms  in  Ae  do  not 
integrate  to  zero  in  this  case.  Thus,  no  advantage  over  the 
one-to-one  case  is  evident. 

The  one-third-to-one  resonance  case  will  eliminate  the 


two-thirds-to-one  resonance  terms  in  Ae  and  vise  versa. 
However,  either  choice  results  in  no  resonance  terms  in  A  i 
or  A/1,  which  contradicts  the  performance  objectives. 

Any  integer  resonance  above  two-to-one  results  in 
nearly  all  changes  being  equal  to  zero  over  one  orbit.  The 
five  terms  in  Ae  which  cannot  be  integrated  to  zero  can  be 
set  to  zero  by  appropriate  choice  o-f  orientation  angles. 
These  resonances,  then,  represent  orbit-keeping  strategies. 
Consequently,  they  are  not  of  interest  for  present 
purposes. 

A  consideration  of  specific  cases  is  provided  in  the 


next  chapter 


V.  Results  and  Discussion 


Resul ts 


For 

the 

single 

resonance  case,  many 

w 

a 

m 

4 

U 

were  tested, 

ranging 

from 

the 

trivialy  simple  to 

the 

arbitrari ly 

complicated. 

These 

cases  were  tested 

by 

numerical ly 

integrating  Eqs  29-32  and  comparing  the  results  with  the 
analytic  solutions,  Eqs  39-42.  The  numerical  integration  was 
accomplished  using  a  predictor — corrector  technique  contained 
in  a  subroutine  Haming,  contained  in  Appendix  B.  The 
Earth-Canonical  system  of  units  was  used  and  a  value  of  D 
was  taken  as  4.65  E-06.  For  the  test  cases,  the  semimajor 

axis  was  taken  to  be  five  DU® . 

Motion  in  the  Ecliptic  Plane. 

The  most  trivial  cases  are  those  in  which  the 
inclination  is  zero.  This  leaves  the  longitude  of  the 
ascending  node  undefined,  and  the  problem  reduces  to  three 
equations. 

A 

lEi.DQi.Q9a.  The  simplest  case  is  the  one  in  which  H 

A 

is  parallel  to  k,  and  the  nutation  angle  is  one-half  pi 
radians.  This  leads  to  the  analytical  expressions: 
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Aa  =  -  a  DTsm(0-r) 


(44) 


Ae  =  A 1.  =  o 


(45) 


Obviously,  th«  maximum  change  in  A  a  is  when 


fl-Y  -  -  trA 


(46) 


This  is  verified  in  Fig  9,  where  the  line  represents  the 

change  in  A  a  over  time  and  the  "x"  indicates  the  point 

predicted  by  Eq  44. 

A 

Tumbling.  Tumbling  is  simply  spinning  with  H 
parallel  to  j  or  i.  The  case  of  H  parallel  to  i  is  trivial, 
since  the  edge  of  the  sail  is  sunward  and  no  thrust  is 
produced.  This  case  was  numerically  verified.  For  the 
non-trivial  case  we  have  the  analytic  equations: 

Aa  =  -3  a3/2DTsin(0-y)  <47) 

i 

Ae  =  O  (48> 

Ai  *  -a/2DTsin(0-Y'-) 


(49) 


The  intsresting  result  is  that  both  Aa  and  Ai  vary 
with  the  sine  of  the  phase  angle  difference.  So  -for  this 
case,  changing  Aa  implies  changing  Ai.  This  case  was 
numerically  verified  and  the  results  are  shown  in  Figs  10 
where  the  line  is  the  variation  over  time  and  "x"  marks  the 
point  predicted  by  Eqs  47-49. 

Coning.  The  simple  caning  case  has  H  directed  at 
the  sun  and  a  nutation  angle  o-f  one-fourth  pi  radians.  The 
resulting  analytical  equations  are: 

Aa  =  -H  a^DTW, <-y)  (50) 
H 

Ae  -  o  (5i> 

Ai  =  -  -FL  a/z  DTsin  (52) 
H 

These  equations  indicate  two  important  results.  First, 
by  selecting  the  phase  angle  difference  to  be  in  the  fourth 
quadrant,  both  As  end  Ai  can  be  positive.  Second,  and 
most  important,  by  selection  of  the  magnitude  of  the  phase 
angle  difference,  it  is  possible  to  trade  off  A  a  for  Ai 
and  vice  versa.  Any  combination  of  A  a  and  Ai  within  the 
magnitudes  given  is  possible.  This  includes,  of  course, 
selecting  only  A  *  or  only  Ai.  Two  cases  were  verified 
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numerically:  the  maximum  Ai  case  and  the  maximum  A  a  case 
both  showed  excellent  agreement  between  analytic  and  numeric 
solutions.  The  results  are  contained  in  Figs  11. 

Motion  in  the  Inclined  Plane. 

The  next  level  of  complexity  is  to  consider  cases  in 
the  inclined  plane.  This  adds  equations  -for  the  change  in 
longitude  of  the  ascending  node.  Initially,  an  inclination 
of  thirty  degrees  was  considered  and  the  simple  cases,  with 
some  variations,  were  considered. 

Spinning.  Two  spinning  cases  were  considered:  case 

A  A  A 

one  with  H  parallel  to  k  and  case  two  with  H  perpendicular 
to  s  and  j.  The  analytical  equations  for  case  one  are: 


Aa  =  -  coszi  aViDTsin(0-r) 


(53) 


Ae  -  AL=  AXX  *  O 


(54) 


As  expected,  case  one  showed  substantial  change  in  A  a 
(.274E-02  DU$)  and  no  other  changes.  These  are  the 
analytical  equations  for  case  two: 


Aq  =  (3  V?  +2)gV:!DTcos(d-H 
6 


(55) 


Ae  -  O  (56) 


AL  *  -_3  a^DTsin  {#-*) 

16 


A  -Cl  -  3a'/;1DT  cos(^-y) 
8 


Cass  two  featured  a  phase  angle  difference  of  negative 


one-half  pi 

radians 

and  so 

showed  significant 

change  in 

A  i 

(. 14E-03) 

and 

no 

other 

changes. 

This 

verifies 

the 

importance 

of 

the 

phase 

difference 

in  determining 

the 

changes  in  the  orbital  parameters. 

Tumbling.  Two  tumbling  cases  were  considered:  case 

A  A  A  A 

one  with  H  parallel  to  j  and  case  two  with  H  parallel  to  i. 
The  analytical  equations  for  case  one  are: 


(57) 


(58) 


Aa^  a3/z DT  3cos(0-yO-5sin(0'V') 

8  i 


(59) 


Ae  =  o 


(60) 


AL  = 


DT  sin 

(6 


(61) 
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V*  A  A  A  . 


<d  _  n  _  »  J  . 
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AH  =  g/1DTl3cos(^-y')-V3’sin(^-y')] 

8 


Case  one,  with  a  phase  angle  difference  of  negative 
one-half  pi  radians,  showed  significant  changes  in  A  a 
( .  228E— 2  DU*>,  Ai  (.137E-3)  and  Ail  (.158E-3).  This 
indicates  a  possibility  that  similar  cases  could  be  used  to 
provide  a  variety  of  changes  in  A  a,  Ai  and  Ail,  depending 
on  mission  requirements.  Case  two  analytical  equations  are: 


Aa  =  <1  aDT  cos  Lp-t) 
16 


A  e*0 


Ai  =  -  3  a2,  DT  sin(^-y-) 

32 


Ail  =  _3_a/zDTcos  (p-v) 
\6 


Case  two  was  tested  with  the  same  phase  angle 
difference  as  the  previous  cases,  providing  significant 
change  in  Ai  (.685E-4)  only.  This  strategy  provides  the 
option  of  changing  i  or  a  and  il  ,  or  combinations  of  all 


three. 


QSQlQSf.  Four  cases  of  coning  were  considered:  two 

A  A 

each  of  case  one  (H  parallel  to  s)  and  two  each  o-f  case  two 

A  A 

(H  parallel  to  i ) .  In  case  one,  phase  angle  differences 
were  changed  to  test  the  case  for  maximizing  A  a  and  Ai, 
and  likewise  for  case  two.  The  analytical  equations  for 
case  one  are: 


Aa  =  VJ  a^DTcos^-y)  (67) 

A  e®o  <68) 

Ai.®  Vi"  a?1  DTsin  <69> 

Ail®  V6  a'/*DTcos(^-?')  (70) 

8 

The  equations  for  case  two  are: 


Aa  =  [876  +  i3Vr)  a/z  DT  cos  {</>-  f) 
6H  J 


(71) 


-  45  - 


A  .V.'i 


A 


(73) 


Ai  -  - 1 5  Vz  a'^DTsinC^-^) 

12.8 

A  SI  =  15  V3C  a‘/zDTcos(8-^)  <74> 

6H- 

Case  two  proved  slightly  superior  -for  changing  Ai, 
(.121E-3  vs  . 112E-3)  with  no  other  orbital  parameters 
changing.  Case  two  also  proved  superior  in  the  second  phase 
angle  difference  case,  giving  superior  performance  in  A  a 
< . 217E— 2  vs  ■  646E— 3  DU.),  and  A-fl  (.242E-3  vs  .224E-3)  with 
no  other  orbital  parameters  changing.  Note  that  maximizing 
a  involves  maximizing  -Cl  .  Trade-offs  can  be  made  between  A  a 
and  Ai.  Agreement  between  numeric  and  analytic  results  was 
excellent. 

Arbitrary  Cases. Three  cases  were  tested  in  which 
phase,  orientation,  inclination  and  nutation  angles  were 
randomly  determined.  Agreement  between  numeric  and  analytic 
results  was  excellent.  One  such  case  is  displayed  in  Figs 
12. 

Nutation  Angle. 


The 

search 

for  the 

optimum  nutation 

angle 

was 

accomplished 

on 

a  case  by 

case 

basis. 

The 

case  of 

A 

H 

pointing 

at 

the 

sun,  inclination 

equals 

zero, 

yields 

an 

optimum  nutation  angle  of  approximately  35.2  degrees.  The 
change  in  A  a  as  a  function  of  @  for  this  case  is  shown  in 

A 

Fig  13.  For  the  30  degree  inclined  case  with  H  parallel  to 

A 

i,  the  optimum  nutation  angle  is  35.2  degrees  for  increasing 
Aa  but  is  40.2  degrees  for  maximizing  A_fL .  For  the  30 

A  A 

degree  inclined  orbit  with  H  parallel  to  s,  the  optimum 
nutation  angle  is  again  35.2  degrees  for  changing  A  a  and 
40.2  degrees  for  maximizing  Ai.  These  are  surprising 
results,  to  say  the  least. 

Limits  of  Validity. 

A  study  of  the  variation  of  A  a  with  a  was  done  to  find 
at  what  point  the  perturbation  assumptions  became  invalid. 
The  study  indicated  a  maximum  change  in  A  a  of  0.25  DUS  at 
-10  DUg  for  the  spinning  case.  This  will  serve  as  an  upper 
bound  on  the  usefulness  of  these  equations. 


Discussion. 

Conclusions. 

The  spinning  case  is  preferred  for  missions  where 
changing  the  semi major  axis  is  the  only  requirement.  Few, 
if  any,  such  missions  exist.  This  implies  that  correct 
inclination  was  achieved  at  insertion  into  orbit,  which 
implies  that  impulsive  power  sources  were  used.  The  cost  of 
such  sources  can  be  avoided  by  using  the  sails  capability 


for  inclination  changes. 


The  tumbling  case  is  to  be  avoided  unless  mission 
requirements  call  for  the  ratio  of  Aa/Ai  inherent  in  such 
a  strategy.  This  is  unlikely. 

The  coning  cases  are  the  only  ones  which  provide  the 
opportunity  to  trade  off  Aa  for  Ai.  Any  ratio  of  Aa/Ai 
desired  can  be  produced  by  adjusting  the  difference  in  the 
phase  angles.  Coning  motion  strategies  thus  compose  a  set 
of  strategies  of  considerable  practical  interest. 

The  3S.2  degree  nutation  angle  results  are  suspect 
until  a  more  detailed  study  can  be  accomplished.  It  is 
possible  that  for  the  differences  in  the  cases  studied  that 
the  differences  in  the  nutation  angles  were  small  enough  to 
escape  notice. 

Directions  for  Further  Study 

More  questions  have  been  raised  than  have  been 
answered.  Besides  the  obvious  point  of  investigating  the 
optimum  nutation  angle  in  more  detail,  a  host  of  other 
topics  suggest  themselves. 

The  many-orbit  case  needs  investigation.  This  involves 
establishing  an  inertial  reference  frame  in  the  plane  of  the 
ecliptic  and  updating  many  of  the  angles  assumed  constants 
in  this  study.  Complex  missions  can  thus  be  analyzed. 
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Including  eccentricities  not  equal  to  zero  will  add 
much  to  the  generality  of  the  results,  but  the  valiant  soul 
who  attempts  this  should  be  forwarned:  the  equations  in 
A  (esinco  )  and  ^  (ecoso;  )  will  be  approximately  -forty 
hand-written  pages  each.  This  author  made  excursions  in 
pursuit  o-f  this  grail  when  the  complexity  of  it  was  still 


cloaked 

in  a 

sophomoric  error. 

Note  that 

this 

level  of 

effort 

only 

gains  enough  terms 

to  make 

consideration  of 

eccentricities  up  to  0.01  possible.  One  might 

be  wel 1 

advi sed 

to 

find  a  dependable 

numerical 

method 

for  such 

pursuits. 

Other  improvements  would  include  adding  the  effects  of 
eclipse  to  the  study.  This  would  consist  of  developing  an 
algorithm  which  would  give  times  of  eclipse,  integrating  up 
to  the  boundary  and  then  coasting  through  the  eclipse. 
Another  improvement  would  be  considering  sails  which  are 
less  than  perfectly  specular. 

Lastly,  there  are  infinitely  many  cases  to  be 
investigated  with  the  equations  provided  herein.  Variation 
of  the  orbital  parameters  with  i  or  with  the  resonance 
number  would  be  a  good  jumping  off  spot  for  the  curious. 

A 

Can  •  function  be  found  which  gives  the  best  angle  between  H 

A 

and  s?  The  possibilities  are  endless. 


Appendix  A 


GgngralPerturbation  Solutions 

For  simplicity  of  notation,  the  -following  two  -functions 
are  defined: 


C  (a,  b)  *  cos  ((an  +  bv)T+ajz$  +b>0  -  cos  (ajrf-b?') 


<A1) 


an  +  by 


S(a,b) 


sinf(an  +  bv)T+a^+b>')  -sin(a^-bf) 


an  +  bf 


(A2) 


$ 


where  a  and  b  are  constants  here. 


Since  these  -functions  typically  occur  in  pairs  such  as 
C(ayb),  C<a,-b)  it  can  easily  be  seen  that  the  resonance 
condition  is  simply: 


a  n 
b 


(A3) 


where  a/b  is  the  resonance  number,  designated  “x"  in  the 
text. 


The  general  perturbation  solutions  follow: 
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Acl=2(x/zDi^\a  Lace.  ,0  -*■  C(  1,-i)  +  C(  1,3)  +  C(t,-3)] 

8 

+A  t  S C  1,-3)  -  SCi,3)  +  S(l,-  0  -  S(l,  lil  +A  [  c(l ,o)*'A(  CC  1,2)  + 

8  2 

C(l,-2))]  +^y[  3S(l,-l)  +  S(l,  l)  ♦S(l,3)-*S(l,-3l]+/£s[c(lJ-0- 
8  8  ^ 

cCiJi)  +  c(iJ-3)-cCi^)]  [5(1,0)  +  yz(s(i,-z)+s(i,2))],) 

, 

-^J1dz(-/d|[s(lr3)  -  S(lJ3)+S^J-l)-S(lll)]+^[-c(lj)-CM 

8  8 

+  c(li3)  +  c(lr3)]-/4[s(lrz)-5(iiz)]+/|4[c(lrl)-c(lj) 

+  C(l,-3)-cO,3)]  +  A[5(lJ-t)  +  S(lj)  -5(  1,-3)“  5(1,3)]  +  A 

8  W 

[c(i,-z) -C(l,2)]/  -^ZJld3(  ^[C ( 1,0)  +  Vi  (  C  ('1, 2)  +  C  ( 1,-Z ) )] 

+/4[5(iJ-z)-S(i,z)]  +^[0(1,1)  +C(l,-  0]  +^[5(1,0)  +  /2(s(iJ-z) 

+  s(i,z))]  ^5[c(i;-z)-c(i,z)]  +|6[ s(i,-0  +  5  0,0]) -£<^3 

(-/^,[S(lrZ)  “50,2)1  ^[-C(l,o)  +  ^(c(l,Z)  +  C(l,-z))]  - 

^CsCiri)-  S(i,0]  +/^v[c(irz)-C(i,2)]+/^’5r 5(i,o)  -  Vz( 

5(l,-2)+  S(l,z))]  ♦/£[c(l/-l)-C(i,l)])  +  [-C(ljl)-C(lj-l) 

z  8 

+  C  ( 1, 3)  +  C0,-3)]  +^[3(5(1,-!)  -5(1,1))  -S(l,-3)+  s(l,3)]  - 

fez  T-C  (l,o)  +  ,/z(ctl,2)  +  c(l,'2))]  +^[s(l,-l)  +  5(l,/)-S(l,-3) 
'I  8 

-5(l,3)]  +  A E3  (c(  1,-1)  -  C(i,  1))  +  C0,3)  -  c(l,-3)]  */§i  [5(1,0)  ~ 
'/i(s(l,-2)  +  S(i,l)}])  +  J3  (  ACc(i,l)  ->c(l,-l)]  +^[s(l,-l)  - 


5(|,|)]  +43C(i,o)  +4i/[s(i,-|)  +  s(ij)] 


(Ai) 


Ac  =  a/xD[  (£(/&[ 3(c(^i )  +  C(2,-i))  +  C (2,3)  +  <2(2, -3)] 

+/4[S(V3)-  SU,3)  +  5(2J-0-S(iJl)3  +4*[c(2,o)  +  Vi( 

'16  ^ 

C(2,2)  +  C(2,-2))]  +  aC^S^J)  +35(0,3)  +  »/i(3{s(z-})  + 

8  / 

5(2,  l)}  +  S(2,-3)  +  S ( 3 ))]  +/^5[~3(c(oi,3)  +  c(o,  0)  +  Vz  ( 

-cM  +  c(v3)  -  C(2,l)  +  C(2,-l))]  +/^[3(T+S(o,z)) 

+  5(2,0)  +  14(5(2/2)  +5(2,2))])  -2<u(-4,[s(v3)-S(2,3) 

5(2, -I )  -  5(2,1)]  +/4[c(2,-3)  +  C(2,s)-C(2,  l)  -c(2,-))] 

16 

-43[s(2r2)  -5(2,2)]  +4<[-3(c(o,3)  +  c(oji))  +  '/i(-C(2,3)  + 

8  8 

C  (2,-3)  - C  (2,  l)  +  C(2,-l)]  ■+  4s[ 3 ^ 5(o, l)  -  S (0,3))  +  '/z[s (2,- ) ) 

+  S(2,l)-S(2,-3)-S(2,3)]  +^[-30(0,2)+  ’/2  ( C  (2,-2 )  -  C  (2, 2 ) )]) 
+  dA-Alc (2,-3)  +C (2,3)- c(2,  l ) -  C (2,- 1 )]  +42[3(  SU-O- 

M6  16 

5(2,1))-  S(2,-3)  +  5(2,3)]  +  a[-C(2,o)  +  '/lfc(2,2)+c(2,-2))] 

H 

+£,[3(5(0,!)-  5(0,3))  +  '/2(5(2,-l) +5(2,1)- 5(2, -3)-5(2,3))] 
+|s[-5C(o,l)  +  3C(o,3)  +  72(3(-c(2,l) +  C(2,-0UC(2,3)- 
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C (2,-3')')!  L  3 ( T -  5(0,2))  +  S(z,o)  -  Vx{  sLz-z)  +  S(z,z))\ ) 

+  Zd1d3(^,[c(2,0)  +  {/x(c(l,z)  +C(zr2))]+/^i[S(zrZ)-5(2J2)] 
+  ^3[C(V')  +  C(2J)]  +^[3(T+S(o,2.))+  Sh.,0)+'A.{s(zrl) 

+  su;2)]+^5[-3c(o,i)  +yi(dvz)-c(^))]  ^[35(0,0  + 
J4(s(irl)  +  s(2,0)])  -2dzd3(^,[S(Z/-2)-Sfo,2.)]  -^[-0(2,0). 

+  !4(c(5L,2)+C(z,-2))]  ts(vi)-  S(z,l)]  +  ^¥  [-3c(oJz)  + 

lA(c(ir2.)-C(2.Jz))]+^5[3(T-  S(o,2.)]  +  S(2.,o)  -  ,/z(S(2J-Z)  + 

S(2.,Z)]  +^[-30(0,0  +  ^C(vO-C(^l})])+d3(/4[c(2rn  + 

H 

C(a.i)]  ^[s(2,-l)-S(Z,l)]  *AC(i,o)  +  A,[3SM  + 

+S(z,l)]  ^5[-3C(o,i)  +  /zfc(z,-l)-C(i,l))]  +^[3T+S(2,o)]1 


(A5) 


Al=  a/z Dl  d,  V /^7L-C(ij)  +  C(t-l)-C(i<3)  +  c(i;-3)]  +^[ 
-C(i,i)-C(l1-l)  +  C(i,3)  +  C(i,-3)])-2j,cl1(/|7[s(i,-i)tS(i,l)- 
5  ( i,  -3)  -  S(  1,3)]  +  d  [c  ( 1,-2)  -  c(  1,2.)]  )  +2  d,  J,  [  c  (\,-z)  -C  (1,2)] 

+^8[s(i,-i)  +  S(i,i)]j  -ZJzJ3(/^7[s(i,o)-  J4  (s(i,-2)  +5(1,2))] 

+  a[c(|,-|)-C(|j|)])+c^(  A[3(c(|,-0-C(I,I))+c(i,3)- 
C(l,-3)]  +  *?B[s(l,o)  -  }i(s(i,-2)  +  S(i,z))])  +  dj  (^[cCl,-!}  - 


An  -  aVlD  CSC  L  L  d'l yS7l  5(  \r3)-S(\5)  +S(lrn-5(lJ/)]  - 

^s[c(  1,0)  +  !/i(c(  u)  +c(l,- -z))])  -2  I )  - 

C(i,-l)  +  C(i,3)  +  C(i,-3)]  +  A,[-c(i,o)  +  !4(c(i,2)-c(i,-z)U 
+  dz  (^[3(5(1,-!)  -S(ljl))  -  S(t,-3)  +5(l,3)]  +ySs[.-C(  1,0) 

+  J4(cd,z)+c(i,-z))]J  +  z«uU[s(ir2)-s(.,d  -/ggl 
c(lj)  +-  c(l,-i)]) -ZeJida^/^rcCi.o)  +  >i(cft,z)  +  cC(,-Z))l 
+/8Cs(i,-,)-S(l,0])  +  A\{ySX  S(y)-S(l,l)]-</S8c(l(o))J 


Appendix  B 
Hami_nq 


-  !*  v  *  **>  «-  pi  f  **  -I  {  n  .<  t,  ) 

'naming  is  an  ordinary  dtferential  equation  integrator, 
x  is  the  independent  variable,  in  our  case:  time, 
y  is  the  state  vector. 

f  is  the  right  hand  sides  of  the  differential  equations 

to  be  integrated,  f  is  explicity  defined  in  a  user-written 
subroutine  called  rhs. 
errest  is  an  error  estimate. 

n  is  the  number  of  ordinary  differential  equations, 
h  is  the  step  size. 


ham  must  be  in  the  calling  program  and  in  rhs. 

The  first  array  parameter  of  y  nnj  f  :n  haai  is  the 

'•'umber  of  DDEs  . 

T h e  fir  3 1  cal!  o n  Haming  i  s  with  •  = 1 .  n • :  t  -  i ,  y  a  t  5  rt 
values  ,  n*--sor  GPEs.  h  is  o ;  oh  cv  .  ond  ' 

call  o  v  i  1  t :  m  so  .joh  •  '  :•  ;  a  on  an  i  ,-<j 
'  R  of  1  ■<  5 


common  /ham/  , 

•col  »  a.. ■/.•/, :.;i 
!  f  i  nxt )  10.0,13, 

13  xo  =  x 

h  h  =  h  /  2 . 3 
call  rhsU) 
do  40  1  =  2,4 
=  x  +  hh 


“  i 


ar  r 


st  i  13  )  .  n  ,  h 


do  2 

0  i  = 

1 

20 

y<  i  . 

1  1  = 

■/( 

call 

rhsi 

1  ) 

::  - 

n  +  h  h 

do  30  i  = 

1 

30 

y(  i . 

i )  = 

y< 

40 

cal  1 

r  h  s  ( 

i ) 

j  sw 

=  -10 

50 

1  sw 

=  1 

do  1 

20  1 

= 

h  h  = 

J’  1  . 

l  > 

i 

>  f( 

;  , 

i  f  ! 

abs  (  ! 

h: 

i  sw 

=  0 

70 

y<  1  , 

2  )  = 

hh 

hh  =  y{  i  ,  1  )  +  h*»  f(!.l'  +  4  ,;df  I  •  ,  2  )  + 
If!  ab  s  (  t  hh -y  (  i  ,  3  )  )  /  ( v  (  :  .  3  ) "' h  )  )  .  1 1 .  to " 
i  sw  =  3 


90 

y(  1  , 

3  ) 

=  hh 

hh  * 

y ! 

i  .  1  5  * 

h* 

!  3, 

.0*7 

t 

+ 

3 .0*f  (  1 

,4 

) 

/  3 

if! 

abs 

1  (  h  h  -  y ! 

i  . 

4 ) )/<?< 

i  sw 

=■  0 

1  10 

y!  1  . 

4  ) 

=  hh 

-  5.0*f!i 
)  go  to  7: 


,  3  >  )  /  2 . 0 
go  to  90 


123  continue 
X  =  xo 

do  1 3JJ  1  =  2,4 
x  =  x  +  h 
130  call  rhs ( 1  ) 

if(isw)  140.140,150 
140  jsw  =  jsw  +  1 

if(jsw)  50.280,200 
150  x  ■  xo 


do  160  1  =  l.n 

160 

er  r est (  i )  =  0.0 

nxt  = 

1 

go  to 

230 

190 

Jaw  = 

2 

n:<t  = 

i  a  b  s  (  n  x  t ) 

200 

X  =  X 

♦  h 

np  1  = 

mod ( nxt , 4 ) 

i 

go  to 

(210.230) . 

!  sw 

210 

go  to 

(270,270,270, 

220) .nxt 

220 

1  sw  = 

2 

230 

nm2  = 

mod ' np 1 , 4  ) 

+ 

1 

nml  = 

nod  (  n;n2 , 4  ) 

f 

1 

npo  = 

mod ( nm 1 , 4 ) 

+ 

1 

do  240  i  =  l.n 

"it,  n,n  2  '  =  v'i.noi)  +  4.0*  h*<  2  .0*1*  <  1  ,  np 

1  *  2  .  J-'f  :  i  .  nn2  >  >  ,  2,0 

240  ■»  v  •.  i  .  r...r2  -  ■}  .  ’  or  r  os  i 

call  rhs i rp  t } 
cc  250  i  *  l . n 

i  ,  no!  ;  =  <  2  ..0**y<  i  ,  npo>  -  :,i\  i  . ->;n2)  •*  2 
1  *  Z  .  0*f  <  i  ,n?o  I  -  '  i  ,  n  ni  1  ;  )  )  /  3  ,  .2 

errsst  (  i  =  f  ;  I  ,  nn2  )  -  y(  )  .  no  1  ) 

250  yM.npl)  =  yil.npi)  -  0.07  42001553  *  err 
jo  to  (  230 , 270  )  ,  j  3W 
250  call  r hs ( np  1  ) 

270  r.xt  =  npl 
230  return 
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pressure  force.  The  equations  of  motion  are  expanded  from  the  Lagrange  planetary 
equations  in  their  acceleration  form  ?nd  are  solved  assuming  small  changes  in  the 
orbital  elements  over  the  period  of  one  orbit. 

Resonances  are  observed  between  the  mean  motion  of  the  sail  and  its 
precession  rate.  The  one-to-one  resonance  case  is  examined  and  numerical  methods 
are  employed  to  verify  the  analytic  results  for  the  circular  orbit  case. 
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